CARMA Processes driven by 
Non-Gaussian Noise 



Robert Stelzer" 



We present an outline of the theory of certain Levy-driven, multivariate stochastic 
processes, where the processes are represented by rational transfer functions (Continuous- 
time AutoRegressive Moving Average or CARMA models) and their applications in 
non-Gaussian time series modelling. We discuss in detail their definition, their spec- 
tral representation, the equivalence to linear state space models and further proper- 
ties like the second order structure and the tail behaviour under a heavy-tailed input. 
Furthermore, we study the estimation of the parameters using quasi-maximum like- 
lihood estimates for the auto-regressive and moving average parameters, as well as 
how to estimate the driving Levy process. 



1 Introduction 

In many applications an observer (scientist, engineer, analyst) is confronted with series of data 
originating from one or more physical variables of interest over time. Thus, he has an observed 
(multivariate) time series and will often either be interested in removing (measurement) noise 
to extract the signal more clearly or in modelling the observed process, including its random 
components. 

In both situations stochastic models may very well be appropriate. This is clear when one is 
mainly interested in removing noise, but when intending to model the observed value it is also 
very often appropriate to enrich a physical model by a random component to capture fluctuations 
and shortcomings of the physical model. The driving stochastic process (the "noise") may have 
interest on its own (as is the case with economic models), or it has to be modelled well to extract 
the interesting information as well as possible (e.g., as is common practice in telecommunication 
links) 
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The easiest way to obtain a model with randomness for the variables of interest would be to 
assume that all observed values are independent and identically distributed (iid) random variables 
or that they follow a physical model plus iid noise. However, in most series observed consecutive 
values are heavily dependent and thus more sophisticated models are needed. A flexible but at 
the same time very tractable class of models is given by linear random processes. In the discrete 
time setting these models are well-known as autoregressive moving average (ARMA) processes 
and they are given in terms of a general order linear difference equation where an iid noisy input 
sequence introduces all randomness. The latter is also referred to as linear filtering of a white 
noise. 

In many situations it is more appropriate to specify a model in continuous time rather than in 
discrete time. These include high-frequency data, irregularly spaced data, missing observations 
or situations when estimation and inference at various frequencies is to be carried out. Moreover, 
many physical models are formulated in continuous time and, hence, such an approach is often 
more natural. 

In the following we consider linear random processes in continuous time, referred to as con- 
tinuous time autoregressive moving average (CARMA) processes. Intuitively, they are given as 
the solution to a higher order system of linear differential equations with a stochastic process as 
the input, which can be seen as linearly filtering the random input. 

One important question is which random input to take in the continuous time set-up. Clearly, 
the random process should correspond in some sense to the idea of white noise. Understanding 
the latter in the strict sense means using independent increments, in the weak sense it means 
uncorrelated increments and so the variance has to be finite. Recall that for random variables un- 
correlatedness is equivalent to independence only if the random variables are Gaussian, i.e. they 
have a normal distribution. A linear random process driven by Gaussian white noise has again 
Gaussian distributions. However, in many situations it is not appropriate to assume Gaussianity 
of the variables of interest, since the observed time series often exhibit features like skewness 
or heavy-tails (i.e. very high or low values are far more likely to occur than in the Gaussian 
setting), which contradict the Gaussian assumption. Demanding uncorrelated but not necessarily 
independent increments does not lead to a nice class of processes nor to nice theoretical results. 

Hence, a good modelling strategy where the resulting process is reasonably tractable and the 
driving process' probability distribution is allowed to have "fat tails" is to demand that the ran- 
dom input shall have independent as well as stationary increments, i.e. increments over time 
intervals of the same length have the same distribution. They then have a time homogeneity fea- 
ture and resemble the iid noise of the discrete time set-up. The resulting class of possible driving 
processes are the so-called Levy processes, which have been studied in detail and form a both 
highly versatile and highly tractable family. An interesting feature is that linear processes driven 
by general Levy processes may exhibit jumps and thus allow the modelling of abrupt changes, 
whereas Gaussian linear processes have continuous sample paths. 

In the remainder of this paper we proceed as follows. First, we introduce Levy processes 
in detail. Thereafter, we give a proper definition of CARMA processes, discuss their relation to 
linear filtering via a stochastic Fourier (spectral) representation and summarize central properties 
of CARMA processes. Next, we briefly explain the equivalence to linear state space models 
and the relation to stochastic control and signal processing. Finally, we discuss the statistical 
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estimation of the parameters and the underlying Levy process and conclude with some additional 
remarks. 

Throughout we will focus on developing the main ideas for CARMA processes. For more 
mathematical details as well as comprehensive references we refer the interested reader to the 
original literature especially the works Q\, &G3l» E-Q, El, 13 ,10] and [Q. For a 
historic perspective the monograph ll47ll may be interesting as well as lll9ll which is the first 
paper where Gaussian CARMA processes appeared under the name of Gaussian processes with 
rational spectral density. 



2 Levy processes 

A Levy process L = (L r ) ?eR + is a stochastic process with independent and stationary increments. 
In the following we consider only Levy processes taking values in the m-dimensional vector 
space W" (with R the real numbers and m some positive integer). Note that a stochastic process 
(Xf )?eR+ can be either seen as a family of random variables indexed by the positive real numbers 
IR + or as a random function mapping the positive real numbers to W 11 . More precisely we have 
the following definition: 

Definition 2.1. An M. m -valued stochastic process L = (L ? ) feR + is called Levy process if 

• Lq = a. s., 

• L tl —L tl ,L t3 —L t2 ,. . . , L tn —L tn _ x are independent for all n G N and t\,t2, ■ ■ ■ ,t n G E+ with 
<h <h < ... <t n , 

• Lf+h — L t = L s+ f l — L s for all s,t,h G IR ("=" denoting equality in distribution), 

p 

• L is continuous in probability, i.e. for all s G M + we have L t — L s —>• as t — > s. 

It can be shown (cf. ll53"ll for a detailed proof) that the class of Levy processes can be character- 
ized fully at the level of "characteristic functions", which we now introduce. Let < -, • > indicate 
the natural inner product in M. m and X is an M m -valued random variable, then its characteristic 
function is defined as \yx(u) = E (e' <u - x> )). The characteristic function of a Levy process can 
always be represented in the Levy -Khintchine form 

E (V<"' L,) ) = exp{t\]/ L (u)}, W > 0, u G R' n , (2.1) 

with 

y L (u) = i(Y,u)--(u,ZGu) + J {j^-l-i(u,x)l m {\\x\\)v{dx), (2.2) 

where 7 G Eg isamxm positive semi-definite matrix and v is a measure on M. m that satisfies 
v({0}) = and / (||x|| 2 A 1) v(dx) < °°. The measure V is referred to as the Levy measure of L 

R m 
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and ||jc|| 2 A 1 is short for min{ ||x|| 2 , 1}. Finally, l^(x) generically denotes the indicator function 
of a set A, i.e. the function which is one if x is an element of A and zero otherwise. Together 
(7, Eg, v) are referred to as the characteristic triplet of L. 

Regarding the paths of a Levy process, i.e. the "curve of L as a function of time t, it can be 
shown that without loss of generality, a Levy process may be assumed to be right continuous and 
have left limits. 

It should be noted that many well-known stochastic processes are Levy processes. Examples 
are Brownian motion, also referred to as the Wiener process or "Gaussian white noise", the 
Poisson process, which has jumps of size one and remains constant in between the jumps, which 
occur after iid exponentially distributed waiting times, and a-stable Levy motions, sometimes 
called Levy flights. Compound Poisson processes are Poisson processes where the fixed jump 
size one is replaced by random iid jump sizes independent of the interarrival times of the jumps. 
It can be shown that all Levy processes arise as limits of such compound Poisson processes. 

A better understanding of what Levy processes really are is provided by the Levy-Ito decom- 
position of their paths. It states that a Levy process is the sum of the deterministic linear function 
yt, a Brownian motion with covariance matrix Lq, the sum of the big jumps which form a com- 
pound Poisson process and the compensated sum of the small jumps (i.e. the sum of the small 
jumps minus their expected value). The quantity v(A) gives for any measurable set A C W n the 
expected number of jumps with size in A occurring in a time interval of length one. In Figure 
Q]a univariate Levy process which is the sum of the linear function t, in this case with 7= 2, a 
standard Brownian motion, with Eg = 1, and a Poisson process, with v({l}) = 1, v(IR\{l}) = 
is depicted together with its individual components. 

Whenever / R m(||.x|| A l)v(dx) < °°, we can replace the compensated sum of small jumps simply 
by the sum of the small jumps adjusting also the slope of the deterministic component. We have 
actually already done this in Figure \T\ where the resulting slope of the deterministic function is 
7— J^xv(dx) = 1. If v(R) < 00, we have finitely many jumps in any bounded time interval and 
the jumps form actually a compound Poisson process. Otherwise, we have infinitely, but count- 
ably many jumps in any bounded time interval. The reason why we have in general a component 
referred to as "the compensated sum of the jumps" (i.e., it results from a certain limiting pro- 
cedure see e.g. [HH]) is that in general the jumps are not summable. This is equivalent to the fact 
that the paths have infinite variation, like Brownian motion. Infinite variation intuitively means 
that the curve described by the stochastic process over finite time intervals has an infinite length. 
Clearly, this means that the fluctuations of the process over small time intervals are rather vivid. 

In Figures [2] and |3] you can see simulations of different pure jump Levy processes, i.e. in these 
cases 7=0 and E = 0. So there is neither a deterministic drift nor a Brownian motion present. 
All these processes have infinite activity, i.e. infinitely many jumps in any time interval. Figure[2] 
depicts a so-called normal inverse Gaussian Levy process which has heavier tails than a Brownian 
motion, but still is rather tame, because it has finite moments of all orders, i.e. E( \L t \ r ) < °° for all 
t , r e M + and also some exponential moments. In contrast to this the stable processes of Figure [3] 
have very heavy tails, because they do not have a finite variance and the 0.5-stable processes does 
not even have a finite mean. Whereas the NIG and 1.5-stable processes have infinite variation, 
the small jumps of the 0.5-stable Levy process are summable. 

Most of the time we will work with Levy processes defined on the whole real line, i.e. in- 
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Figure 1: A Levy process and its components: The complete Levy process is depicted in the 
lower right display. In the upper row the deterministic drift component is depicted on 
the left and the standard Brownian motion component on the right. The left display in 
the middle row shows the standard (rate one) Poisson component and the right one the 
Brownian motion and the deterministic component added together. In the last row on 
the left the Brownian component plus the Poisson jumps are depicted. 
Note that the scaling of the y-axis is different in the individual plots. 
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Figure 2: Simulation of a Normal Inverse Gaussian (NIG) Levy process. 
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1.5-Stable Process 
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Figure 3: Simulations of stable Levy processes. A 1.5-stable Levy process is depicted in the 
upper row and a 0.5 -stable in the lower one. 



7 



dexed by R not R + . They are obtained by taking two independent copies of a Levy process and 
reflecting one copy at the origin. 

For detailed expositions on Levy processes we refer to |Q]], |0], [37] or [53]. 



3 Definition of CARMA processes and spectral 
representation 

On the intuitive level one wants to be able to interpret a <i-dimensional CARMA(p,q) process Y 
as the stationary solution to the p-th order linear differential equation 

P{D)Y, = (DP +A l D?- 1 + ... +A p )Y t (3.1) 
= (BoDV + BiDP- 1 + . . . + B q )DL t = Q(D)DL t , (3.2) 

where the driving input L is an m-dimensional Levy process, D denotes differentiation with 
respect to t, and the coefficients A i, . .. ,A p are d x d matrices and Bq, . . . ,B q are d x m matrices. 
The polynomials P(z) = z p +A i z p ^ 1 +... +A p and Q(z) = Bq^+Biz 9 ' 1 +...+B q with z G 
C are referred to as the auto-regressive and moving average polynomial, respectively. Finally, 
p,q G N are the auto-regressive and moving average order. 

However, the paths of non-deterministic Levy processes are not differentiable and so the above 
equation cannot directly provide a rigorous mathematical definition. Let us briefly consider the 
case (p,q) = (1,0) in which case the resulting process is actually called an Ornstein-Uhlenbeck 
(OU) process. In the univariate case it is given by the differential equation 

DY t =aY t +DL t 

where a is a real number. So what we basically want is that the change of Y over an infinitesimal 
time interval is a times the current value of the process times the "length of the infinitesimal time 
interval" plus the change of the Levy process over the infinitesimal time interval. Rephrasing this 
idea in the precise language of stochastic differential equations (see e.g. I48I1 ) we obtain 

dY t = aY t dt + dL t . 

Using the theory of stochastic differential equations (SDEs) it is easy to see that this SDE has a 
unique solution given by 

Y, = e at Y + e at [ e- as dL s . 
Jo 

For general orders (p, q) one could to some extent use a similar reasoning to arrive at a precise 
definition of CARMA processes. However, we shall take a more elegant route. First note that 
the differential operators on the auto-regressive side of (13.11) act like integration operators on the 
moving average side. Hence, they offset the differential operators of the moving average side 
acting on the Levy process. Since Levy processes are not differentiable, we effectively have to 
integrate at least as often as we differentiate to be able to make sense of (13.11) . Hence, a necessary 
condition ensuring the proper existence of CARMA processes is p > q. 
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In order to obtain a rigorous definition of CARMA processes our strategy here shall be to 
switch from the time domain to the frequency domain where the main tool is the following 
spectral representation of a Levy process. Here and in the following we denote by A* for a 
matrix (or vector) A the Hermitian, i.e. the complex conjugate transposed matrix. 



Theorem 3.1 (042Q). Let (L t ) te ^ be a square integrable m-dimensional Levy process with mean 
E[L\\ = (which implies E[L t ] = for all t) and variance E[L\L*] = at t = 1. Then there 
exists a unique m-dimensional random orthogonal measure with spectral measure such 



that E [4>l(A)] = Ofor any bounded Borel set A, Fi(dt) = ^dt and 

L t = / — & L (dn),teR. 

J-oo 111 



The random measure is uniquely determined by 



e -i\ia _ e -i\ib 



*L([a,b))= I — dL^ (3.3) 



for all — oo < a < b < 



The random orthogonal measure <J»^ can intuitively be thought of as the "Fourier transform" 
of the Levy process. If L t is a Brownian motion, then [0,f)) is again a Brownian motion. For 
general Levy processes rather little can be said about the properties of 4>l. For example, it is 
known that <I>l has second-order stationary and uncorrelated increments, but the increments are 
neither independent nor stationary in a strict sense, see 02711 . 

In the spectral domain we can now interpret differentiation (and integration) as linear filtering 
noting that a formal interchange of differentiation and integration gives "DL t = e'^' 4>L(J/i)". 
It can be shown that the resulting process is well-defined whenever the linear filter is square 
integrable. Thus we obtain as definition for "7(f) = P(D)- l Q(D)DL(t)": 



Definition 3.2 (CARMA Process, 114211 ). Let L = (L t ) te ^ be a two-sided square integrable m- 
dimensional Levy-process with E[L\\ = and E[L\L*] = E^. A d-dimensional Levy-driven con- 
tinuous time autoregressive moving average process (Y t ) te ^_ of order (p,q) with p,q e No and 
p> q (CARMA(p, q) process) is defined as 

oo 

Y t = J e ifit P(in)~ l Q(in)& L (dn), teR, where (3.4) 

— oo 

P(z): = I mZ p +A 1 z p - 1 + ...+A p , 
Q(z): = B z q + B lZ q - 1 + ....+ B q and 

<i>L is the Levy orthogonal random measure of Theorem I3.il Here Aj 6 M m (R), j = 1, and 
Bj e M d)tn (R) are matrices satisfying B q ^ and^V(P) :={zeC: det(P(z)) = 0} C K\{0} + M 
(i.e. the autoregressive polynomial has no zeros on the complex axis). 
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Referring to the explicit construction of the random orthogonal measure 4>^, one can easily 
show that the above defined CARMA processes are necessarily stationary (in the strict sense, 
i.e. the distributions are left unchanged by a time shift). Since by construction any CARMA 
process in the sense of Definition 13.21 has a finite variance, it is also weakly stationary, i.e. the 
second-order moment structure (the variance and autocovariances) are left unchanged by time 
shifts. 

Although the definition of CARMA processes via a spectral representation is elegant and help- 
ful in many theoretical considerations, it is not really usable in applications, as alone simulating 
a CARMA process from this representation would be a tedious and problematic task. However, 
luckily we have the following result. 



Theorem 3.3 (State Space Representation, 114211 ). Let the Levy process L and P, Q be as before. 
Define the following coefficient matrices: 

p-j-i 

• Pp-j = ~ I Aifip-j-i + B q _j, j = 0, 1 , . . . , q, J8i = . . . = P p - q - 1 = 



7=1 



JS* = (JS 1 *,J8 2 *,...,JS;) and A 







-Ap 



-A p -\ 



Denote by G t = (G\ t , . . . , G* f )* a pd -dimensional process and assume that ^V(P) :={zGC: 
det(P(z)) = 0} C (— °°, 0) + iR - the open right half of the complex plane. Then 

dG t =AG t dt + PdL t (3.5) 

has a unique stationary solution G given by 

G t = f e A{ - t ~ s ^dL s , teR. (3.6) 

It holds that 

/oo 
e l ^P(iii)- l Q{i^ L {dil) =Y t ,te K. 
-oo 

So the first d-components of G are the CARMA process Y. 

A CARMA process satisfying Jf{P) := {z G C : det(P(z)) = 0} c (-°o,0) + /R is called 
causal, because as shown above the value at a time t only depends on the Levy process up to time 
t, it is a function of (L s ) se (_ 00 t y In other words a causal CARMA process is fully determined 
by values in the past. Whenever the condition jV(P) := {z G C : det(P(z)) = 0} C (-°°, 0) + iR 
is not satisfied, Y t also depends on future values of the Levy process. In many applications, 
where it is clear that all we see today can only be influenced by what happened up to now, 
one only considers causal processes as appropriate models. However there are also applications 
where non-causal processes are useful. For example, if we want to stochastically model the water 
level in a river and think of t as describing the location along the river, both the water levels 
downstream (in the "future") and upstream (in the "past") may influence the water level at a 
certain point. Note that in this paper we only discuss stationary CARMA processes. In some 
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applications (e.g. control) it is often adequate to consider non-stationary (non-stable) systems. 
Then the roots det(P(z)) = 0} in the set (— °°,0) +iM. describe the stable and causal part of the 
system and the remaining roots describe the non-stable part. 

Theorem 13.31 allows us to treat a causal CARMA process as a solution to the stochastic differ- 
ential equation (13.51) and thus we can apply all the available results for SDEs. In particular, tasks 
like simulation of a causal CARMA process are straightforward and easily implemented. How- 
ever, the above result allows us also to get rid of another restriction. So far we could only define 
CARMA processes driven by Levy processes with finite second moments and thus we could so 
far not have e.g. CARMA processes driven by a- stable Levy processes. However, general the- 
ory on multidimensional Ornstein-Uhlenbeck processes (see 06h and [Hi]]) tells us that (13.61) is 



the unique stationary solution to (13.51) as soon as the Levy process has only a finite logarithmic 
moment. 



Definition 3.4 (Causal CARMA Process, [42]). Let L = (L ? ) ?e R be an m- dimensional Levy pro- 
cess satisfying 



J In ||*|| v(dx) <oo, (3.7) 

p,q G No with q < p, and further A\,A2, . . . -,A pi G Mj(WL), Bq,B\, . . .,B q G M^ m (M), where Bq ^ 
0. Define the matrices A, /3 and the polynomial P as in Theorem \3.3\ and assume o(A) = C 
(—oo,0) + iM. Then the d-dimensional process 

Y t = (I d ,0,...,0)G t (3.8) 

where G t = f t _ 00 e A ( t ~ s >fidL s is the unique stationary solution to dG t = AG t dt + f5dL t is called 
causal CARMA(p,g) process. 

G is referred to as the state space representation. 

A natural question is clearly whether one can also extend the definition of CARMA processes 
via the spectral representation to the case with infinite variance. For so-called regularly varying 
Levy processes with finite mean and thus especially for a-stable Levy processes with cc G (1,2) 



a result like Theorem 13. II has been established in [27]. However, the non-finite variance case is 
distinctly different, as a limit of integrals has to be taken and the random orthogonal measure is 
replaced by an object which is - strictly speaking - not even a measure anymore. In that paper 
a definition of CARMA processes with regularly varying Levy input analogous to Definition 
13 .21 has been given and it has been shown that the resulting processes coincide with the causal 
CARMA processes when both definitions apply. Observe that processes with infinite variance 
are not only of academic interest, but that they have important applications, for instance, in 



network data modelling (cf. 114311 and [50, 51]). In 112911 CARMA processes driven by a-stable 



Levy processes have been successfully used to model electricity prices. 

4 Properties 

In this section we explain and summarise various properties of (causal) CARMA processes. 
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4.1 Second Order Structure 



Recall that for convenience we have assumed that the driving Levy process and thus the CARMA 
process has mean zero. Looking at the "defining" differential equations, it is clear that if E[L\ ) = 
/l then the CARMA process is defined as the one driven by L\ — \it plus A~ l B q /j, which is then 
the mean of the CARMA process. 



Proposition 4.1 ([42]). Let Y be a (causal) CARMA process driven by a Levy process L with 



finite second moments and set E^ = var(Li). 

1. The CARMA process Y has autocovariance function: 



f e l>xn i , 

cov(Y l+h j t ) = j —p(ivy l Qm^LQ{W(.pmydv, 



with h G 



2. IfY is a causal CARMA process, its state space representation G has the following second 
order structure: 



var(G f ) = J e Au pL L p*e A * u du 



o 



Avar(G f )+var(G,)A* = -/3E L /3* 

cov{G t+h ,G t ) = e Ah var(G t ),h>0. 



Since we are only considering stationary CARMA processes, the moments above do not de- 
pend on t. 

Since Y is given by the first d components of G the second order structure of G implies im- 
mediately alternative formulae for the second order structure of Y . In particular, it shows that the 
autocovariance function always decays like a matrix exponential for h — > °°. 



4.2 Distribution 

Another nice feature is that in principle the distribution of a CARMA process at fixed times as 
well as the higher dimensional marginal distributions, e.g. the joint distribution of the process at 
two (or n) different points in time, is explicitly known in terms of the characteristic function. The 
reason is that all these distributions are infinitely divisible and that their Levy-Khintchine triplet 
is known in terms of the Levy-Khintchine triplet of the driving Levy process. We state this in 
detail for the stationary distribution in the causal case. 

Proposition 4.2 ([H^]). IfL has characteristic triplet (7, E, v), then the stationary distribution of 
the state space representation G of a causal CARMA process is infinitely divisible with charac- 
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teristic triplet (}^,Lg, Vq), where 

oo 

•7g = [e As p 7 ds + J J ^*i3x[l[o I i](||^*J8x||) - l [0 ,i](IWI)] v(dr)^, 

K m 

f' OO 

• Eg= / e As pLp*e A * s ds, 

Jo 

• V £(fl) = f° / l s (e A ' 5 /3x)v(^)^ 

JO >/R m 

for all Borel sets B C R^ rf . 
/« of/?er words 

£(>' G '>) = exp j /(y- M )_I( a ,E^)+ | (^-1-/( M ^)1 [0)1] (|W|)Vg(^) L (4.1) 
for all u e R pd . 

Projection onto the first d coordinates gives the characteristic triplet of the stationary distri- 
bution of Y. It should, however, be noted that typically the distribution of the CARMA process 
does not belong to any special family of distributions even if one starts with especially nice Levy 
processes. 



4.3 Dependence Structure 

An important property of multivariate stochastic processes, is how their future evolution depends 
on the past. Suppose that one stands at a given point in time and one disposes of sufficient data at 
that point to determine the evolution from that point on, also given knowledge of the input from 
that point onwards. A Markov process is a stochastic process, for which the future only depends 
on the current value and not anymore on the past values (all their information is subsumed in the 
current value). For a Markov process it - so to speak - only matters where we are now not were 
we came from. If this characterising property does not only hold at all fixed times, but also at 
certain random times called stopping times, we speak of a strong Markov process. 

Proposition 4.3 (F42I1). The state space representation G of a causal CARMA process is a strong 
Markov process. 

Intuitively it is desirable in many applications that the farther away observations are in time, the 
less dependent they should be. Usually, one even wants that very far away observations should be 
basically independent. This idea is mathematically formalized in various concepts of asymptotic 
independence often referred to as some form of "mixing". 

A comparably weak result which, however, applies to any CARMA process is the following. 

Proposition 4.4 ([28]). Any stationary CARMA process is mixing. 
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Mixing implies ergodicity, i.e. empirically determined moments from the time series converge 
to the true moments if more and more data is collected. So time averages converge to ensemble 
averages. This is very important for statistical estimation of CARMA processes, as it implies 
typically that estimators are consistent (i.e. the estimators converge to the correct value when 
more and more data is collected). 

Typically, one also wants to know the errors of estimators which can be derived from distri- 
butional limit results like asymptotic normality. To obtain such results a stronger more uniform 
notion of asymptotic independence is needed, which is called strong mixing. Typically, one can 
best establish it for a Markov process. 

Proposition 4.5 ( 114211 ). For a causal CARMA process withE(\\L\ \\ r ) < 00 for some r>0the state 
space representation G and the CARMA process Y are strongly mixing, both with exponentially 
decaying mixing coefficients. 



4.4 Sample Path Properties 

Next we look at the sample path properties of a CARMA process. 
Proposition 4.6 ([42]). 

• The sample paths of a CARMA(p, q) process Y with p > q + 1 are (p — q—\ )-times differ- 
entiable and for a causal CARMA process it holds that 

d' 

—Y t = G i+htl i=l,2,...,p-q-l. 

• If p = q + 1 and the driving Levy process has a non-zero Levy measure V satisfying 
v(Bq (M. d \{0})) 7^ 0, then the paths of a CARMA process exhibit jumps and the jumps 
sizes are given by AY t :=Y t — Y t = BoAL t . 

• If the driving Levy process L is a Brownian motion, then the sample paths ofY are con- 
tinuous and (p — q—l )-times continuously dijferentiable, provided p > q+\. 

For examples of the paths of CARMA processes driven by an NIG Levy process see Figure HI 



4.5 Tail Behaviour 

As already stated in the introduction one may want to move away from Gaussian models, be- 
cause extreme (i.e. very low and/or high) observations are far more likely than in a Gaussian 
distribution. One says that the tails (of the distribution) are heavier than Gaussian ones. Very 
often it appears also reasonable to use models which are "heavy-tailed" in the sense that only a 
limited number of moments exists, i.e. Z?(||X|| r ) exists only for low values of r. Mathematically 
it is then convenient to use the concept of regular variation (see 1E2I1 or [49, \5^\ for compre- 
hensive introductions in relation to extreme value theory). Roughly speaking this means that the 
tails behave like a power function when one is far from the centre of the distribution. A random 
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10 20 30 40 50 60 70 80 90 100 

Time 

Figure 4: A CARMA(1,0) process driven by an NIG Levy process having discontinuous paths 
is shown in the upper display and a CARMA(2,0) process driven by the same Levy 
process having continuous paths in the lower one. 
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variable X is regularly varying if P( \\X\\ > x) behaves comparably to x~ a for some a > and 
big values of x. In Il44n (see also [[23(1 in the univariate case) it is shown that under a very mild 
non-degeneracy condition a CARMA process driven by a regularly varying Levy process is again 
regularly varying with the same index a. Hence, it is straightforward to construct heavy-tailed 
CARMA processes when applications call for such features. 

In the univariate case the tail behaviour of CARMA processes is also understood in certain 
non-Gaussian situations, where one has lighter tails than regularly varying ones (see 1241 12511'). 



5 State space models 

We have defined the causal CARMA process using a so-called state space representation and 
we have noted that the state space representation G is made up of the CARMA process Y and 
its derivatives as long as they exist. Hence, causal CARMA processes may be viewed as special 
state space models driven by Levy processes. In fact, any state space model can also be realized 
as a CARMA process, as will be shown now. 

We start with a precise definition of state space models. 

Definition 5.1. Let L be an m-dimensional Levy process and 

A G Mn(M), BeM N>m {R), CeM d ^ N (R). 

A general (N ,d) -dimensional continuous time state space model driven by L with parameters 
A,5,C is a solution of 

the state equation dX t =AX t dt + Bdh t 
and the observation equation Y t =CX t . 

X is called the state process and Y the output process. 

Note that the state process is Af-dimensional whereas the output process is J-dimensional. 
Sufficient conditions for the existence of a unique causal stationary solution of the state equa- 
tion are given by indicates the "real part" of a complex number or function) 

9t(A v ) < 0, A v , v = 1, . . . ,N, being the eigenvalues of A 

and L having finite second moments. 

It can easily be shown by integration that X satisfies 

Xt = e A{ts) Xs + J' e A{t-u) Bdu _ 

Likewise, the stationary output process Y satisfies 

Y t = f Ce A ^BdU, 

J — CO 
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Its spectral density, the Fourier transform of the autocovariance function, is given by 
f Y (co) = -Lc(ico-A) r l BZ L B T (-ia-A T ) r l C 7 '. 

From Definition 13 .41 it is obvious that a CARMA process is a (p d, d) -dimensional state space 
model driven by an m-dimensional Levy process. The following theorem states that also the 
converse is true. 



Theorem 5.2 ( 15511 '). The stationary solution Y of the multivariate state space model (A,5,C, L) 
is an h-driven CARMA process with autoregressive polynomial P and moving average polyno- 
mial Q if and only if 

C(zI N -A)- 1 B = P(z)- l Q(z), VzGC. 
For any (A, 5, C) there exist P, Q such that the above equation is satisfied and vice versa. 

In reality we typically do not observe some variables of interest continuously, but only at a 
discrete set of points in time. Let us assume that we sample the process at an equidistant time 

(h) 

grid with grid length h > and denote by Y„ := Y n f, for neZ the sampled observations of a 
state space process. 
It is easy to see that 



X? } = ***X?> + I"* e A ^BdL U: (5.2) 

J(n-l)h 



Yi h) =CX { n h) (5.1) 

(■nh 

!(n-l)h 

which immediately shows that Y„ is the output process of a discrete time (N : d) -dimensional 

state space model driven by the Af-dimensional iid noise ( J/^n^ e A ^ nh ^ u " > Bdh u ) 

It is well-known that any (N,d) -dimensional state space model in discrete time is an ARMA 
process. Combining this with Theorem 15 .21 tells us that any equidistantly sampled CARMA pro- 
cess is an ARMA process. This observation will be the basis for estimating CARMA para- 
meters in the next section, where we will need a considerable refinement of this result. 

In many applications the sampling frequency is quite high, i.e. h is very small. Thus it is 



important to understand how F^ behaves as h — > which has been investigated in fl 1 8M . 

As we only observe the process Y in a state space model, an important question is what can 
be said about the state process X based on the observations. Hence, we want to reconstruct or 
"estimate" the latent process X as good as possible. This procedure is also referred to as filtering. 



For Gaussian state space models the easily implementable Kalman filter (see e.g. 111111 ) is optimal 
both from a variance point as well as a distributional point of view. For non-Gaussian state space 
models with finite variance the very same procedure, now typically called linear filtering, gives 
an "estimate" of the latent process which is the linear (in the observations) "estimate" with the 
lowest variance. However, it is typically not the "estimate" with the minimal variance and not a 
conditional expectation. Thus, there are more involved filtering techniques, like particle filtering 
(see e.g. I120ll ). which are better. 
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State space models, mainly Gaussian ones, are also heavily used in stochastic control (see 113011 
and references therein for a comprehensive overview) and signal processing (see [38,[39|,|45|], for 
instance). In both areas one is sometimes dealing with data for which a Gaussianity assumption 
is not really appropriate due to skewedness, excess kurtosis or heavy-tailedness. Clearly, in such 
situations Levy-driven state space models or equivalently CARMA processes should be appeal- 
ing. Going into the details of the usage in control is beyond the scope of this paper, but it seems 
worthwhile to mention that there are two uses of state space models in control. Sometimes one 
assumes that one has some random input which is then "controlled" by the state space model, 
so the the state space model acts as the controller. In contrast to this sometimes the output of 
the state space model is regarded as the natural output of some system on which an additional 
controller is acting to ensure that the output meets certain requirements. 



6 Statistical Estimation 

In this section we discuss ways to estimate the parameters of a CARMA process and its driving 
Levy process. First we address the estimation of the autoregressive and moving average para- 
meters. Due to parametrisation issues explained later on, we formally do this for Levy-driven 
continuous time state space models, as defined in the previous section. In the univariate case 
quasi-maximum likelihood estimation of CARMA processes is comprehensively studied in [ITvh . 



6.1 Quasi-maximum likelihood estimation 

We assume that we observe the process Y at discrete, equally spaced times 

Y { n h) :=Y nh , «GZ, h>0. 
Furthermore, we define the linear innovations £™ by 

Jh)_ v (h)_ p V (A) 
c n — x n r n— 1 x n : 

where P n \ denotes the orthogonal projection onto span : — °° < V < n \, i.e. the linear 
space spanned by the observations until time (n — \)h. From the construction it is immediate that 

(h) 

(e„ ) ne z is a white noise sequence, i.e. it has mean zero, a constant variance and is uncorrelated. 
The construction implies that one can only sensibly speak of linear innovations when the driving 
Levy process has finite second moments. Thus we will demand the latter for the remainder of 
this section. 

Theorem 6.1 (JHH]). Assume the eigenvalues X\ , . . . , Xn of the matrix A are pairwise distinct and 
define complex numbers 4>i, 4>2, . . . , <&n by 



l-4> lZ -4> 2 z 2 -...-4>^= Y\ 



N 

' WzeC. 



v=l 



1 - e- Xvh z 
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Then there exist 0i , 02, . . . , ®n-\ in M c j(C) such that 



v (*) *,vW d> -p(*) j. j. 4_<a 

*n -^l*,,.! -5, + 0i£ J? _ 1 + ... + 0j V -i£ )1 _;v+i 

holds. 

Hence, w a weafc ARMA (Af,iV — l) process. 

This result suggests that one could estimate simply the ARMA coefficients of the sampled 
process and then transfer these estimates to estimates of the CARMA coefficients. However, to 
estimate a CARMA process it is not sufficient to estimate an ARMA process, because not all 
ARMA processes can be embedded in a CARMA process. There are ARMA processes which 
cannot arise as equidistantly sampled CARMA processes. The way out is carry out the "ARMA 
estimation" in the CARMA parameter space. 

Since we are going to use a quasi-maximum likelihood approach and have discretely sampled 
observations, all possible models considered in the estimation have to be distinguishable based 
only on the second-order properties of the sampled process. 

Definition 6.2 (Identifiability). A collection of continuous time stochastic processes (Y#, # G 0) 
is identifiable if for any #i^^2 the two processes Y#j andY$ 2 have different spectral densities. 

It is h-identifiable, h > 0, if for any ^ ^ the two processes Y^ and Y^ have different 
spectral densities. 

We assume that our parametrisation is given by a compact parameter space 0cM ? with some 
q G N and a mapping 

Here, A$ is the N x N matrix of our Definition 15 . 1 1 dependent on the parameters -& and likewise 
for B$,C$ and L#. 

We need to ensure that our parametrisation is minimal regarding the dimensions, since a fixed 
output process can result from artificially arbitrarily high-dimensional state space models. 

Assumption PI (Minimality). For all # G the triple (A^.B^.C^) is minimal in the sense that 
if 

C(zI m -A)- l B = C$(zI N -A$)- l B i ) 

then m>N must be true. 

Assumption P2 (Eigenvalues). For all # G the eigenvalues of A$ are pairwise distinct and 
contained in the strip 

{zeC:-7t/h<3(z)<7t/h}. 

We want to use a parametrisation for the continuous time state space model, but need to ensure 
that it is /z-identifiable. The following theorem provides easy-to-check criteria. 



Theorem 6.3 ( Il56ll ). Assume that the parametrisation \j/ : D & i— > (A^,5^,C^,L,j) is 
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• identifiable, 



• minimal 

• and satisfies the eigenvalue condition. 

Then the corresponding collection of output processes {Y#, i?6 0} is h-identifiable. 

The quasi-maximum likelihood (QML) estimator is now obtained by pretending the observa- 
tions were Gaussian, taking the corresponding likelihood and maximising it. More precisely the 
QML of ■& based on L observations y L = (y i , . . . , yi) (of a CARMA process with parameter t9n) 
is 

=argmax t>e0 «^ (y L ), 
where Jzf^ is the Gaussian likelihood function which is proportional to 

-1/2 



EI det eX P l - \ £ e l,nV#fa,n 



i I 2 , 

i«=l / I n=l 



with 



{Y^ v =y v :l<v<n} 



V» n =E 



r 



Y^ = y v :l<V<n 



So „ are the linear innovations under the model given by # and V$ : „ are their variances or 
the one-step prediction errors. Note that y L are in contrast to this observations of the CARMA 
process with the unknown parameter #o which we are about to estimate. 

Computing the QML estimator is now a straightforward task utilising the Kalman recursions 
and numerically maximising the likelihood. However, since we have not used the true likelihood, 
it is not clear whether the resulting estimators are really sensible in the sense that they converge 
to the true parameters. Luckily, one can show that the estimators are well-behaved. 



Theorem 6.4 (Strong consistency, IdoI ). Assume the parametrisation \]f is continuous. For every 

A L 

sampling interval h > 0, the QML estimator # is strongly consistent, i.e. 

a L 

■& —7- #o o,.s. as L — > oo, 
provided the parametrisation is h-identifiable. 

However, so far we cannot assess the quality of our estimators by confidence intervals etc., 
which is made possible by the following result. 
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Theorem 6.5 (Asymptotic normality, n56\ ). Assume that the driving Levy process satisfies E\ |L^, ( 1 ) 1 1 4+5 < 
°° for some 8 > and that the parametrisation \j/ is three times continuously differentiable. For 
every sampling interval h > 0, the QML estimator ■& is asymptotically normally distributed, i.e. 



\fl\b -^ )4^K(o,n), a = 7(^ ) _1 /(W(^o) 



7 ^ = / im 7^^ ln ^( yL )' 
/(^) = limivar^-ln^(y L ), 

provided the parametrisation is h-identifiable. 

To obtain identifiable parametrisations one uses like in the discrete time case (see ll32ll or ||40ll . 
for instance) so called canonical parametrisations like the echelon state space form. For more 
details on this we refer to [56] . Since such parametrisations are typically available for state space 
models rather than CARMA processes, one normally estimates state space models rather than 
the equivalent CARMA processes. 

Let us finally look at one simulation study. 

A d-dimensional normal inverse Gaussian (NIG) Levy process L (see e.g. |0, 0, S^]) with 
parameters 

8 >0,k>0,/3 GR rf ,AGM+(M) 
is given by a normal mean-variance mixture, i.e. 

Li =i u + VA/3+V 1/2 Ar, 

where N is <i-dimensionally normally distributed with mean zero and variance A and independent 
of 

V ~IG(5/fc,5 2 ) 

which follows a so-called inverse Gaussian distribution dHH]). 

We consider now a bivariate NIG-driven CARMA process with zero mean given by the state 
space form 



dX t = 



01 02 

1 

3 #4 05 

1 
1 



X t dt + 



01 02 
06 07 
#3 + 0506 04 + 0507 

08 09 

09 010 



dL t , 



The parameters are #1, #2, . . . , #10 and the parametrisation is in one of the canonical identifiable 
forms. 
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Figure 5: One realisation of a bivariate NIG-driven CARMA process (upper two displays) and 
the effect of sampling (lower two displays). The linearly interpolated process over the 
time interval [600, 650] resulting from sampling at integer times is shown as the thicker 
line, whereas the thinner line is the true CARMA process. 
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parameter 


sample mean 


sample bias 


sample 
standard deviation 


estimated 
standard deviation 


0i 


-1.0001 


0.0001 


0.0354 


0.0381 


02 


-2.0078 


0.0078 


0.0479 


0.0539 


#3 


1 .005 1 


-0.005 1 


0.1276 


0.1321 


#4 


-2.0068 


0.0068 


0.1009 


0.1202 


#5 


-2.9988 


-0.0012 


0.1587 


0.1820 


#6 


1.0255 


-0.0255 


0.1285 


0.1382 


07 


2.0023 


-0.0023 


0.0987 


0.1061 


08 


0.4723 


-0.0028 


0.0457 


0.0517 


09 


-0.1654 


0.0032 


0.0306 


0.0346 


010 


0.3732 


0.0024 


0.0286 


0.0378 



Table 1: Summary of the results of the simulation study on the QML estimation of a bivari- 
ate NIG-driven CARMA process. The second column states the mean of estimators 
obtained over 350 simulated paths, the third column the resulting bias and the fourth 
column the standard deviation of the obtained estimators. Finally, the last column states 
the standard deviation for the estimators as predicted by the asymptotic normality result 
Theorem [63] 

A simulated path is shown in Figure [5J 

We calculated the QML estimates for this bivariate NIG-driven CARMA process based on 
observations over the time horizon [0, 2000] at integer times and repeated this for 350 different 
simulated paths. The estimation results are summarised in Table [U It shows that the sample 
bias of the obtained estimators in the simulation study is very small and that the sample standard 
deviation is close to the standard deviation predicted by the asymptotic normality result Theorem 
16.51 Actually, the sample standard deviation is always smaller which is nice, as it implies that the 
standard deviation predicted by the asymptotic normality result Theorem 16.51 is a conservative 
estimate. 

6.2 Statistical inference for the driving Levy process 

The above quasi-maximum likelihood approach only allows to estimate the autoregressive and 
moving average parameters as well as the variance of the driving Levy process. However, typ- 
ically we want to estimate many more parameters of the driving Levy process or even first need 
to get an idea to which family the driving Levy process may belong to. To this end one can re- 
construct from the CARMA process the driving Levy process. Typically, the CARMA process is 
only observed at a discrete set of times and then the best we can do is to get approximations of the 
increments of the Levy process. One can then treat the approximate increments as if they were 
the true ones of the Levy process. "Looking" at them one should be able to choose appropriate 
parametric families. By using the approximate increments, as one would use the true ones, in 
maximum likelihood or method of moment based estimation procedures one can do parametric 
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inference for the Levy process. The construction of the approximate increments and their use 



in estimation procedures has been studied in detail in 111411 where it is in particular shown that 
the estimators are good in the sense that they are consistent and asymptotically normal under 
reasonable assumptions when taking appropriate limits. 



It should be noted that the idea to reconstruct the Levy process can already be found in [10] 



or II 1711 . In the following we illustrate this approach for a univariate Ornstein-Uhlenbeck, i.e. a 
CARMA(1,0), process based on [fl6TI and Bill from which all examples and plots are taken. 
Recall that an Ornstein-Uhlenbeck (OU) process is the unique strictly stationary solution to 

dY t = aY t dt + dL t . (6.1) 

where (L r ) feR is a Levy process with £'(ln(max(|Li|, 1))) < °° and autoregressive parameter 
a < 0. The solution of the stochastic differential equation is given explicitly by 

t 

Y t = e a ^Y s + le^dLu. (6.2) 



If the OU process is observed continuously on [0, T], then the integrated form of (16.11) imme- 
diately gives 

Lt = Y t -Y -a I Y s ds. 



o 



(h) 

The increments of the driving Levy process AL„ on the intervals ({n — l)h,nh] with n G N can 
be represented as 



(h) f" h 

'■= L n h-L( n -i)h = Y nh -Y( n _^ h -a Y u du. (6.3) 

J (n— l)h 

What we want, is to approximately reconstruct the sequence AL^ of increments over intervals 
of length h from observations of the CARMA process made over a finer equidistant grid. To this 
end one simply approximates the integral f" n _^ h Y u du by some numerical integration scheme 

needing only the values of the process on this finer grid. Since the approximations of AL^ 
become thus closer and closer to the true increment as the numerical integration scheme becomes 
more exact, Ill4ll derive their asymptotic results when both the observation interval as well as the 
observation frequency goes to infinity. Note that in practice one does not know a so one has to 
estimate it first, which could e.g. be done by the already described quasi-maximum likelihood 
approach. 

Turning to an example, let us consider the OU process given by 

dX t = -0.6X t dt + dL t . (6.4) 
with L being a standardised Gamma process, i.e. L t has density 

y 1 / 2 ^ , 1/2 
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gampdf(x,y,1/sqrt(y)) 
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Figure 6: Probability density of the increments of the standardised Levy process with 7 = 2 and 
the histogram of the estimated increments from one path of the OU process, obtained 
by sampling the process with grid length 0.01. (Source: Bill .') 

and the parameter 7 being set to 2. 

In Bill 100 paths of this OU process on the time interval [0, 5000] have been simulated and 
then the Levy increments over time intervals of unit length have been approximated by sampling 
the OU process over a grid of size h. 

In Figure [6]the histogram of the Levy increments distribution from one path with h = 0.01 is 
shown, together with the true probability density of L\ . 

If one further averages over all one hundred paths which is equivalent to looking at one path 
over a one hundred times longer time horizon, the fit of the histogram to the true density becomes 
visually almost perfect, see Figure |7J 

Based on the approximate Levy increments one can now estimate the parameter 7 by max- 
imum likelihood. Table [2] shows summary statistics of the resulting estimator for different samp- 
ling grid sizes h. The data in the table is based on estimating 7 separately for each of the 100 
simulated paths. 

To conclude, the simulation study illustrates that the recovery of the background driving Levy 
process and the parametric estimation based on the approximate increments works quite well. 
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Figure 7: Probability density of the increments of the standardised Levy process with y = 2 and 
the histogram of the estimated increments for all 100 paths of the OU process, obtained 
by sampling the process with grid length 0.01. (Source: 113 1IU 



Table 2: Estimated parameters of the standardised driving Levy process based on 100 paths on 
[0, 5000] of the Gamma-driven OU process. 



h 


Parameter 


Sample mean 


Sample standard 






of estimator 


deviation of estimator 


0.01 


7 


2.0039 


0.0314 


0.1 


7 


2.0043 


0.0340 


1 


7 


1.9967 


0.0539 
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7 Concluding Remarks 



Finally, we would like to mention that there are other stochastic models like the so-called ECO- 
GARCH process of [33] and B4I1 where CARMA processes are an important ingredient as 
well as extensions of CARMA processes. One extension are fractionally integrated CARMA 
(FICARMA) processes (see [13] and 114 IIP . While CARMA processes have an exponentially 
decaying autocovariance function and thus have always short memory, FICARMA processes ex- 
hibit polynomially decaying autocovariance functions and are thus able to model long memory 
phenomena (see 112 ill or Ii52ll for detailed introductions into the topic of long range dependence). 
However, the paths of FICARMA processes are continuous. A class of processes with possible 
long memory, jumps in the paths and related to CARMA processes are the supOU processes, 
see HLH] and ll26ll . As noted in multivariate supOU processes can be straightforwardly ex- 
tended to obtain so-called supCARMA processes. Long memory is (believed to be) encountered 
in data from many different areas, e.g. finance or telecommunication. Since it is an asymptotic 
property and similar effects in the autocorrelation function might be caused by structural breaks 
(non-stationarity), it is often hardly debated whether there truly is long memory in a time series. 
The first scientific study considering long range dependence properties was looking at the water 
level of the river Nile (see Ii57ll ). 

From the overview on CARMA processes presented in this paper it should not only be clear 
that they are useful in many applications, but also that there are still many questions to be ad- 
dressed in future research. These include alternative estimators to the ones presented here, es- 
timators which work in the heavy-tailed case when one does not have a finite variance or order 
selection, i.e. a theory how to choose the orders (p,q) of the autoregressive and moving average 
polynomial when one fits CARMA processes to observed time series. 
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